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Fluorescent materials play a prominent 
role in the qualitative and quantitative 
measurement of scientific phenomena 
of importance in biotechnology and 
biomedical applications. Photodegradation 
of fluorophores is a process that 
determines the accuracy and sensitivity 
of such measurements.This is the 
motivation for developing methods for 
accurately measuring fluorophore 
photodegradation rates. Recently, 
illumination consisting of short pulses 
has been used to examine the decay 
of photochemical reaction products. 
However, the time resolved measurements 
are difficult to interpret since the 
photodegradation process usually involves 
multiple time scales. The frequency 
domain measurement technique discussed 
here looks at the frequency response 
of a fluorescent sample to a frequency 
modulated illuminating light. The 
photodegradation rate is obtained by 
interpreting the frequency domain 
measurements in terms of traditional 
impedance concepts. In the measurements 
described in this paper, a focused laser 
beam is used to illuminate a sample of 
slowly flowing fluorescent solution. The 
laser beam is assumed to have a Gaussian 
power distribution hence illumination is 
spatially non-uniform in the region of 



interest. The photochemical reaction rates 
depend on power, so they will also vary 
with the position in the beam. However 
in the case of photodegradation of 
fluorophores, the measurement of the 
resulting decrease in fluorescence is given 
in terms of the radiation emitted from 
the entire illuminated region. In this work 
we present a mathematical description 
of the time evolution of the fluorescence 
response integrated over a non-uniformly 
illuminated domain. As a result of our 
analysis, an experimentally accessible 
and tractable mathematical model 
Eq. (19) and Eq. (30) is obtained from a 
more fundamental description given by 
Eq. (4) and Eq. (5). The model is used to 
create a functional form for fitting 
experimental measurements from a lock-in 
amplifier. 
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1. Introduction 

The photodegradation of fluorophores is an impor- 
tant process that determines the accuracy of many bio- 
logical assays [1], the efficacy of photodynamic thera- 
py [2], and the phototoxicity of drugs. The measure- 
ment of the apparent rate of photodegradation was tra- 
ditionally performed by measuring the decrease in 
fluorescence emission under constant illumination [3]. 
Recently illumination consisting of short pulses has 



been used to examine the decay of photochemical reac- 
tion products [4]. The time resolved measurements are 
difficult to interpret since the photo-degradation 
process usually involves multiple time scales. The fre- 
quency domain measurement technique looks at the 
response for each frequency of modulation of the illu- 
minating light. Harmonic equations are simpler to solve 
and, thus, interpretation of frequency domain measure- 
ments can be performed in terms of traditional imped- 
ance concepts [5]. In the measurements described 
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below, a focused laser beam is used to illuminate the 
sample. The power distribution of a Gaussian laser 
beam is strongly dependent on spatial location; 
therefore, the rates of photochemical reactions, which 
depend on power, will also depend on location in the 
beam. Since absorption of radiation leads to local heat- 
ing and concentration gradients, convective and 
diffusive mass transport will be present in the vicinity 
of the focused beam. A slow flow is imposed on the 
fluorophore solution in order to dominate and thus min- 
imize the effect due to uncontrollable convective and 
diffusive mass transport. In the case of photodegrada- 
tion of fluorophores, measurement of the resulting 
decrease in fluorescence usually measures the radiation 
emitted from the entire illuminated region. Therefore 
an analysis method has to be developed which 
describes the time evolution of the integrated fluores- 
cence response from a nonuniformly illuminated 
region. In previous work [6] a frequency domain meas- 
urement technique was interpreted with a model which 
assumed a uniform beam profile. This work extends the 
model to a laser beam with a Gaussian profile. In it we 
derive an experimentally accessible mathematical 
model Eq. (19) and Eq. (30) from first principles Eqs. 
(4) and (5). A dramatic simplification of these equa- 
tions is a consequence of the strongly disparate time 
scales in the kinetics imposed by the slow rate of degra- 
dation relative to the relaxation and absorption rates 
The first section contains a motivation and presentation 
of the equations that are taken as a fundamental 
description of the experimental setting (also see Fig. 1). 
The equations are based on a two-state model with 
excitation and relaxation. In Sec. 2, a perturbation 
scheme employing the disparate time scales is used to 
reduce the system of two partial differential equations, 
describing the kinetic process to a single partial differ- 
ential equation Eq. (14) for the total fluorophore popu- 
lation (to leading order in the perturbation parameter). 
The biggest reduction in mathematical complexity 
comes from writing the total measured fluorescence in 
terms of the total fluorophore population averaged over 
the Gaussian light beam distribution. This formula 
is introduced in Sec. 3, Eq. (9) and there follows the 
derivation of a single ordinary equation for the time evo- 
lution of the spatial average Eq. (30) that completes the 
model. The mathematical arguments that justify the 
manipulations performed in the main body of the paper 
can be found in the appendices. Section 4 discusses the 
time independent case and a functional form for fitting 
experimental measurements from a lock-in amplifier. 
The form, obtained from fitting values predicted from 
the model, shows good agreement with actual data. In a 



forthcoming paper, a harmonic approximation of the 
solution of the model equations will be developed. It 
will be possible then to express the fit parameters in 
terms of the physical parameters in the experiment. 

2. Kinetic Equations 

As a starting point, consider a model for reactions at 
a fixed point in space. This simplified model will be 
used to motivate the details of the analysis. Assuming 
no mass transfer from surrounding fluid. Fig. 1 shows a 
two state model where absorption occurs from the 
ground state with rate k^ , and the excited state popula- 
tion of size TVi, relaxes back to the ground state with a 
rate k^ . Here N^ is the population of the ground state 
molecules. Intact fluorophores change into non-fluor- 
escent species with a rate k^ . The kinetics are summa- 
rized by a set of equations 

^ = -^A+^A (1) 



dt 

dN^ 
dt 



kN-kN-kN^ 



To interpret ^^, Eq. (1) has to be supplemented by a 
model of the photochemical processes that are respon- 
sible for the conversion of the fluorophores into non- 
fluorescent species [6]. In the present manuscript, there 
will be no attempt to interpret the rate constant k^ • The 
discussion below will focus on two problems inherent 



N, 




kr 



Nn 



Fig. 1. A graphical representation of tiie components included in the 
kinetic model discussed in the text. The rate constant k^, represents 
photodegradation that leads to non-fluorescent species. The rate con- 
stants k^ and A:^ represent absorption and relaxation to the ground state 
respectively. Nq and Vj represent the populations of the ground and 
excited states. 
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in using Eq. (1) to describe real experimental condi- 
tions. The first problem is that the laser beam has a 
spatial intensity distribution leading to spatial variation 
of the absorption rate. The second problem is that con- 
vection and diffusion currents will change the local 
concentration of molecules in an uncontrolled manner. 
The second problem will be reduced significantly by 
introducing a flow that will change the concentration of 
fluorophores in a predictable way. 

The absorption rate depends on the power density of 
the incident light. The absorption rate will be written as 
k^ = a J [7], where (7a(cm^) is the molecular absorption 
cross section, and / is the incident photon flux (1/s 
cm^). During the experiment, the laser irradiance 
P(W/cm^) is measured. The irradiance P can be con- 
verted into a photon flux, /, by dividing P by the ener- 
gy per photon. Explicitly / = P^P where the conversion 

factor is given by i^ = /h ' ^^^^ ^ ^^ ^^^ wavelength, 

n is the index of refraction, h is Planck's constant, and 
c is the speed of light in vacuum. This yields the result 
k^ = o^P^P = k^P, The constant k^ relates the laser 
irradiance to the molecular absorption rate. 

Figure 2 shows the geometric relation of the laser 
beam, the cuvette holding the sample, the flow, and the 
detector. The laser beam is incident along the z axis 
which points out of the plane of the paper. The irradi- 
ance of the modulated beam that illuminates the sample 
will be written as: 



P{x, y, t) = Pf^ (x, y) + AP(x, y) cos(cot) . 



(2) 



Poix^y) is the time independent component, and AP(x,y) 
is the amplitude of the modulated component, O) is the 



CO 



radial modulation frequency, — is the frequency in Hz, 

2k 

and / is time. The beam will be assumed to have a 
Gaussian profile [8]. 



fi^.y)- 



1 



-exp( ^— ) 

w 



KW 

P,ix,y) = PJix,y) 
APix.yj) = PJix,y)cosicot) 



(3) 



where Pq is the total power (watts) of the laser beam 
and w is the width of the beam. The function /(x,;^) is 
defined such that its integral over all the x,y plane is 
normalized to 1 . The time independent component and 
the modulation amplitude have the same spatial depend- 




cuvGtts wall 



Fig. 2. The geometric arrangement of cuvette walls (thick lines), 
laser beam propagating in the z direction (out of the plane of the 
paper), the fluid velocity in the x direction, and the detector located 
on the >> axis. The shaded circle represents the Gaussian profile of the 
laser beam. 

ence given by f(x,y). Clearly in this situation, the ab- 
sorption rate will have a strong spatial dependence. 
Consequently reaction products may have concentra- 
tion gradients and uneven heating may cause thermal 
mass transfer. The kinetic model of photodegradation, 
given by Eq. (1), will be extended to include a spatially 
dependent absorption rate and a constant flow of the 
solution containing the reaction product. It will be 
assumed that the flow dominates all mass transport 
(diffusive and convective) into and out of the illuminat- 
ing region and provides a well defined initial concen- 
tration of fluorophore in the illuminated region. 

To account for the spatial variation of absorption and 
the flow of fluid carrying fluorophores past the laser 
light, Eq. (1) is modified to. 



dN^^^dN^ = -K(x,y,t)N,(x,y,t) 
dt dx 

-\-k^N^(x,y,t) 

BN^^^BN^ = k^(x, y, t)N,(x, y, t) 
dt dx 

-k,N,(x,y,t)-kjNf(x,y,t) 



(4) 
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Here k^ {x,y,t) = kcP{x,y,t). The velocity of the fluid, 
denoted by v is assumed to be a constant and in the 
direction of the x axis in Fig 2. Note in contrast to 
Eq. (1) the absorption rate is now a function of space 
and time. To completely specify the solution of these 
equations, the number of particles flowing into the 
entrance of the apparatus must be specified and the 
number of particles at the beginning of the observation 
period (/ = 0) must also be given. 

N,{-L,y,zJ) = N\ N,i-L,y,z,t)=0 for t>0 
N^(x,y,z,0) = N\ N,(x, y, z,0) =0 for x>-L 

(5) 



Equation (5) states that prior to entering the laser beam 
there is no change in the populations of the fluorophore 
states. The distance L is some length along the x axis 
where the laser beam intensity is effectively zero. 
Eq. (5) also gives the initial conditions at / = when the 
laser beam is turned on. 

3. Perturbation Analysis of Kinetic 
Equations 

The faster processes, absorption and relaxation 
dominate the time variation at least initially. Returning 
to the model for the fixed position in Eq. (1), this sug- 
gests that after an initial transient the quasi-equilibrium 
condition holds. 






Substitution of this condition into Eq. (1) results in 
a single equation for the total fluorophore population 

N=Nq + N,. 

dN 
dt 



-KK 



K+K 



■N. 



In this section this reduction is extended to the case 
of non-uniform illumination and is made rigorous. Here 
the coupled pair of partial differential equations, 
Eq. (4), will reduce to a single partial differential equa- 
tion Eq. (14), accurate to first order in a perturbation 
expansion. Justification of the operations leading to 
these simplifications can be found in Appendix A. 

The absorption and the excited state relaxation 
rates are usually of the order 10^ s~\ while the value 
of ^^ is less than 10^ s~\ The constants in Eq. (4) can 
therefore be re-scaled as k^ = k^/ £, k, = kj £ where 
£=10"^ and k^ , k, are of the order of 1. Equation (4) 
can then be rewritten as, 



dN^{x,y,t) ^ ^ dN^{x,y,t) \ 
dt dx J 

-k^{x,y,t)N^{x,y,t)+k^N^{x,y,t) 

dN,(x,yj) ^ jN,(x,yj) y 
dt dx 



(6) 



k^(x,yj)N^(x,y,t)-k^N^(x,y,t)-ek^N^(x,y,t) 

We will treat NQ(x,yt,) and Ni(x,yj) as functions 
of £. Let NQ(x,yJ)=NQ(x,yJ,£)\^=Q and N]^(x,yt) = 
N^ (x,yt,£) |e=obe the leading terms for Nq (x,yt,£) and 
N^ (x,y t, £) respectively in a perturbation expansion for 
£>0 and large time t. 



N^ (x, J, /, £) = N^ (x, J, t) + £N^,ix, y, t) + 0{ e") 
N,ix,yj,e) = N,ix,y,t)+eN,,ix,y,t)+Oie^). 



(7) 



The validity of the expansion in £is discussed further 
in Appendix A. Substituting Eq. (7) into Eq. (6) and 
equating coefficients of equal powers of £ on both sides 
of the equation gives for £ = 



= -k^(x,y,t)N,(x,y,t)+k^N,(x,y,t) 
= k^ix,y,t)N,ix,y,t)-k^N,ix,y,t), 

This implies that as £ ^ 



(8) 



(9) 



Thus the relation between the populations of ground 
and excited states holds in this approximate sense for 
all larger times. (This is not valid for / ~ 0. When the 
illumination is turned on, there is an initial transient 
period where the populations are changing rapidly and 
therefore where the left hand side of Eq. (6) is not pro- 
portionate to £.) Equating the coefficients of the first 
power of ogives: 



dt dx ' " 

-k^(x,y,t)N^^(x,y,t) 

dN,(x,y,t) dN,(x,y,t) y , , - , 
-kMx,y,0-k,Ni(x,y,0- 



(10) 
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Letting N(x,yJ)=No(x,yJ)+ Ni(x,yj), be the total 
fluorophore population in the zero order approxima- 
tion, we can add the two parts of Eq. (10) and obtain 



dN{x,yJ) dNix,y,t) 

+ V- 



dt 



dx 



-k^N^ix,y,t). 



(11) 



The right side of Eq. (11) depends on the zeroth 
order approximation to the population of excited state. 
Using Eq. (9) and the definition of the total fluorophore 
population it is possible to obtain the relation 



N,ix,y,t)= ^°f'/''^ Nix,y,t) 



(12) 



which is correct to zero order in e. Substituting Eq. (12) 
into the right hand side of Eq. (11) gives the following 
equation for A'^(x,>',/) valid up to first order in e 



dN(x,y,t) dN(x,y,t) 



dt 



dx 






N(x,y,t), 



(13) 



Equation (10) has been rewritten as a kinetic equation 
for the total fluorophore number valid to the first order 
in £. Equation (13) can be rewritten further by setting 
b = kj k,, r} = k^b. Eq. (13) becomes 



dNix.yj) _ -^^(-^'-V'O,^. A_v_^^(-^'3^'^) 



dt 



l + bPix,y,t) 



dx 



(14) 



The boundary condition for N{x,yt) at x = -L is 
inherited from those imposed on Nq and A^i at x = - Z in 
Eq. (5) 



N{rL,y,t)=N\ 



(15) 



The initial condition at / = is N{x,y,Qi) =N^ since no 
reaction takes place prior to turning on the illumination. 

The last term in Eq. (14) gives the flow related 
flux of fluorophores into the illuminating region as a 
function of position. The flow of fluorophores into the 
volume illuminated by the laser beam is assumed to be 
uniform over the entire volume. The photodegradation 
depends on the intensity of incident light; consequently 
the concentration of intact fluorophore will not be 
uniform over the illuminated volume. 



4. Spatial Averaging and tlie 
Fluorescent Signal 

The measured fluorescence signal per unit distance 
in the z direction is given by 



F{t) = Ak,^,jJN,ix,y,t)dxdy (16) 



where A represents the characteristics of the optical 
measurement instrument, k^ is the rate of radiative 
decay of the excited state population, and Ni(x,yt) is 
the population of the excited state. (The radiative decay 
contributes to the overall relaxation rate given by k^). 
The population of the excited state can be approximat- 
ed by Eq. (12) leading to 



Fit) = AK^ \\±-^l2D-N{x, y, t)dxdy. 
J J i + bP(x,y,t 



bP{x,y,t 



(17) 



The symbol b was defined previously. Equation (17) 
together with Eq. (14) constitute a description of the 
measurement valid to first order in e. Equation (17) 
suggests that instead of finding the solution N{x,yt) at 
every spatial point it may be better to find the variation 
of the averaged concentration defined by Eq. (18) 



{Nit))^jjfix,y)Nix,y,t)dxdy, (18) 



The average is performed with a weighting factor 
equal to the spatial distribution of the laser irradiance. 
The weighting factor provides a measure of the contri- 
bution of fluorophores at each point in space to the 
total fluorescence signal. It will be shown below that 
Eq. (14) can be converted into an equation for (A^(/)) as 
defined by Eq. (18), and that the fluorescence signal 
can be written as 

F(0 = Ak,,,b[p(t){N(t))-ba(I\i)f{N(t))] (19) 



where to a very good approximation a is a constant and 
P(t) = Po + PiCOs((jD t) = Po + A(t). For actual computa- 
tions, the constant a is evaluated in terms of the time 
independent solution of Eq. 14. In this section, an ordi- 
nary differential equation for the time evolution of 
(A^(/)) will be derived (Eq. (30)). The fluorescent signal 
can therefore be obtained (with small error) by solving 
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Eq. (30) and then using Eq. (19). Together these equa- 
tions represent an experimentally accessible mathemat- 
ical model of the fluorescent signal under excitation by 
a focused laser beam. 

To find an equation for {N(t)}, first multiply Eq. (14) 
by 1 + bP{x,y,t) and obtain (below r represents the 
coordinates x,y). 



(\ + bP{rj))^^^^ = -riP{rj)N{rj) 



dt 
-v(l+Z)P(r,0) 



dNjrj) 
dx 



(20) 



Multiply Eq. (20) by/(x,>') and integrate over the beam 
cross section region, R. Each term will be treated in 
detail. To simplify the notation integrals in the x-y plane 
will be represented by a single integral in r = (x,y) with 
dr = dxdy. The left side of Eq. (20) becomes 

J/(r)(l+^P(r,0)^^rfr 

R 

which can be expanded to 

R R 

Inserting the explicit form of the power function we get 

R 

R 

d(N(t)) d f 

= ^^ + K^o+A(0)-J/(r)MnO/(r)c/r. 

R 

The second part of the above equation will be rewrit- 
ten by introducing the quantity 

\f(r)f(r)N(r,t)dr 
a(t) = ^ . (21) 



jf(r)N(r,t)dr- 



R 



The final result for the left side of Eq. (20) after multi- 
plication and integration is therefore 



^^^ + biP, +Ait))j\ait)lfir)Nir,t)dA (22) 

(Note that the definition given in Eq. (21) will be used 
in deriving Eq. (19)). 

The first term on the right side of Eq. (20) after the 
operations becomes 

-7ljPirj)NirJ)fir)dr 

R 

= -r]{P, + A(0)J/(r)iV(r, t)fir)dr (23) 

R 

= -77a(f)(P„+A(0)(iV(0) 

where a is defined in Eq. (21). The second term on the 
right side of Eq. (20) becomes, after multiplication and 
integration, 

-yj(l + bP(r,ty)^^^f(r)dr 



dx 



dN{r,t) 



=-^\'-^md, 



dx 



-vb{P, + A(0)J/(r)^^^ /(r) Jr. 



To treat the second term of the above equation, we 
introduce a second function defined by 



;^VA - R 



jf(r)m^-^dr 



a'(0 = 



dx 

R 



(24) 



The second term of the right hand side of the trans- 
formed Eq. (20) now reads 

-v(l + b(P, + A(O)a'(0) //(/-) ^^^ dr. (25) 



In Appendix B we will show that the functions (l(t) 
and (Z\t) are, up to a small error, independent oft and 
can be replaced by constants a and a ' respectively. 
The error introduced into equations Eq. (21), Eq. (23), 
and Eq. (25) by this approximation is discussed in 
Appendix C. 
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Collecting terms, Eq. (20) can be written after multi- 
plication by f(x,y) and integration, up to small varying 
functions of t, as 



i\ + baiP,+Ait))) 



d{Nit)) 



dt 



---riaiP, + Ait)){Nit)) 

-y(l + ba\P,+A(t)))l^)(t) 



(26) 



where 



dN_ 
dx 



(0 = j/W 



dN(r, t) 
dx 



■dr. 



(27) 



The final step in the derivation of an equation for 
(A^(/)) consists of finding an approximation to the 
weighted average spatial derivative in Eq. (27). For 
small widths w the major contribution to the integral 
comes from a region in the vicinity of the beam which 
we define byO<|x|,|j|<K*w where we may choose 
the constant k, 0<k<3. If we replace the partial 
derivative in the integrand by a backwards difference 
quotient with step h=- kw the value of A'^ at the left 
endpoint of the interval is A'^ (to a good approxima- 
tion). Thus for w small enough we have. 



dN\^^_{Nit))-N' 



dx 



(0 



KW 



+ 0(1) 



where o(l)is a term that tends to zero as w ^ 0. To 
guarantee an error that is smaller than the dominant 
terms of Eq. (26) we have to require that w is sufficient- 
ly small and K*is large. In fact experimental conditions 
govern the actual choice of w. The following heuristic 
argument governed the choice of K 

We observed that in the vicinity of the origin (x = 0, 
J = ) the spatial derivative of the time independent 
solution, N(x,y) discussed below, has a functional form 
which is almost identical to the beam distribution func- 
tion /(x.j-) as shown by the two overlapping curves in 
Fig. 4. The difference between -f(x,y) and dN(x,y)/dx 
(both normalized by dividing by the maximum of their 
respective absolute values) is magnified 50 times and 
shown in Fig. 4 by the upper curve. The weighted aver- 
age of the un-normalized beam distribution function 
71 ■ width^f{x,y) is equal to 0.5. Therefore the weighted 
average of the spatial derivative will be approximately 
equal to 0.5 times the value at the origin 



(0 = \f{x,y)^;;^{x,y,t)dxdy 



dx 



dx 



\dN_ 

2 dx 



(0,0,0- 

(28) 



The derivative at the origin can be approximated in 
the usual way by the difference of the function N{x,y, t) 
evaluated at two points spanning the origin and divided 
by the distance between the two points. However for 
the purpose of this analysis a more useful approxima- 
tion uses the difference {N{t))-N^ divided by an 
effective distance which reproduces the value of the 
derivative. For the case of the Gaussian beam, the best 
estimate of the derivative at the origin was found 
(using the time independent solution) to be {{N {t)) 
- N^)lw'^0,631, This yields the following approximate 
expression for the average spatial derivative given in 
Eq. (28). 






1.274 w 



(29) 



Using Eq. (29), the approximate kinetic equation for 
(A'^ (/)) becomes 

A(Ar(o) = ^^^(^±^(^(Ar(0) 
dt^ ' \ + ba{P^+A{t)y ' 

_^__(l + .a-(i>+A(0)) 

1.274W (l+to(Po+^(0)) 

(30) 
where A{t) = P^ cos(6>/). 

5, Validation of tlie Model Equations and 
Their Derivation: 

This section discusses the time independent form of 
Eq. (30) (see Eq. (31). Since an analytical form of 
the solution exists for the time independent case of 
Eq. (14), direct comparisons of the fluorescent signal 
obtained from the time independent version of Eq. (30) 
are possible. In both cases the fluorescence signal is 
obtained from Eq. (19) (using the steady state solutions 
of Eq. (14) and Eq. (30). A second important reason for 
considering the time independent solution is that the 
spatial average in Eq. (18) is dominated by the spatial 
average of the time independent solution itself As a 
consequence, the time varying functions a (t) and a\t) 
can be replaced by constants calculated from the time 
independent solution (see B 9 in Appendix B). 

The time independent form of Eq. (14) corresponds 
to a condition of constant illumination, constant flow, 
and constant initial fluorophore concentration. 



197 



Volume 112, Number 4, July- August 2007 

Journal of Research of the National Institute of Standards and Technology 



r\PJ{x,y) 
l + bPJ(x,y) 



N(x,y) = -y 



dx 



(31) 



{N) 



Steady state 



Figure 2 shows the geometrical arrangement of the 
flow channel, the laser beam, and the detector. The 
laser beam is along the z direction, the detector along 
the J axis, and the flow takes place along the x axis. The 
channel dimensions are d^= 0.4 cm along the y axis, 
d ^ = 1 cm along the z axis, and 5 cm along the x axis. 
The velocity of the fluid in the center of the channel is 
set to 0.025 cm/s. No uncertainties are given since 
these values will be used to model the response of the 
system. The width of the Gaussian laser beam is set to 
w = 50 ■ 10^ cm. The laser power is Pq = 0.03 W, the 
vacuum wavelength is A = 488 • 10"^ cm. The number 
of photons per Joule, P^, is calculated by Pc = X n/hc = 
3.314 • 10^^ J~\ where n = 1.35 is the index of refrac- 
tion of water and he is the product of Planck constant 
and the speed of light. The extinction coefficient is 
set to £ = 82000 L/cm mol, yielding a cross section 
cr=1000f;/0.4343A^^ = 3.135- lO-'^cm" where A^^ is 
the Avogadro's number. The average absorption rate is 
given by k^ = oP.PqNk width = 3.97 • 10^ s"\ and the 
relaxation rate is given by the inverse of the measured 
lifetime of the excited state, k^= 1/t=3.33 • 10^ s"\ 
The photodegradation rate k^ is set to 170 s~\ Using 
these values it is possible to calculate the parameters 
that enter into Eq. 3 1, Z) = aPJK, = 3. 1 17 • 10^ cmVJ 
and?7 = ^,Z> = 5.3- lO^^cmVJ. 

The solution of Eq. 3 1 is written as 



1 + - 



WPq 



■N", 



^.26w(l+^«'^o) 



(33) 



The two values of the fluorescent signal (FS) agree 
to within 1 % up to 0. 1 W incident power. The approx- 
imate fluorescence signal given by Eq. (30) and 
Eq. (33) is almost indistinguishable from the result 
given by Eq. (32). It is expected that the time depend- 
ent solutions will also be described correctly by the 
approximate Eq. (30) which together with Eq. (19) 
will be taken as the fundamental description of the 
photodegradation process occurring in the focused 
laser beam. Finally the steady state solution was used to 
find the effective distance to use in the estimate of the 
average spatial derivative given by Eq. (29). The exact 
value of the average spatial derivative and the 
estimated value using Eq. (29) are - 3.932 cm""^ and 
- 3 .906 cm"^ respectively. Thus Eq. (29) should provide 
a reasonable estimate of the average spatial 
derivative even in the time dependent case. 

The formal time dependent solution of Eq. (30) 
can be written as ( see for example 
http://tutorial.math.lamar.edu/terms.aspx) 



(NiO) = N''i^>''"'0 + j; P''"'G(^)dt') (34) 



where 



N(x,y) = Ae'^'''''^ 



dx' 



(32) 



G(0 = 



riaPjt) V \ + ba'Pit) 

\ + baP{t) \21Aw\+baP{t) 

V \ + baP{t) 
\21Aw\ + baP{t) ' 



(35) 



where L = 0.02. Equations (21) and (24) can be used to 
calculate a =6324 and a' = 8462, respectively, using 
Eq. (32). Note that ^ is a constant which is set to 1. 
Figure 3 shows the resulting N(x,y) for the parameters 
given above. The solution given by Eq. (32) can be 
used to predict the time independent fluorescence sig- 
nal by substituting it into Eq. (19). The result can be 
compared to fluorescence calculated from the steady 
state of the approximate kinetic Eq. (30), i.e.. 



For the case of low power, the approximations 
baP{t) « 1, ba'P{t) « 1, and a laser beam power given 
by P(t) =Pq + Pi cos {(O t) yield an approximate solu- 
tion given by Eq. (36). 



(Ar(^)) = 7V^e-t^«^o-/i.274w).^- , 



^^sin(c.O 



(1 + - 



i.274i 



t 



^(7)a/'„+v/1.274w)/'^ ft, 



, 52^sin««') , (36) 

dt). 
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Fig. 3. (a) The steady state solution of Eq. (3 1) given by Eq. (32). The physical parameters that enter into 
the equation are discussed in the text. The function N(x,y) is normalized to 1 . The laser beam is incident 
along the normal to the x-y plane. As expected, the number of intact fluorophores decreases with increas- 
ing laser beam intensity. The decreased concentration of intact fluorophores persists along the x axis due 
to the imposed flow, (b) The steady state jc-derivative of the function N(x,y). The derivative is finite along 
the profile of the laser beam. 
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Fig. 4. The correspondence between the Gaussian beam profile function -fix,y\ shown by the soUd line, and the steady state function 



dN(x,y) 
dx ' 



shown by the dotted line. Both functions are normalized by their maximum values in order to emphasize functional similarities. The difference 
between the two functions, magnified 50 times, is shown by the curve with two valleys. The graphs in Fig. 4a are for>' value of cm, and the 
graphs in Fig. 4b are for>' value of 0.005 cm. 
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Figure 5 shows the time dependent solution at three 
values of the laser power for a modulation frequency of 
2 Hz. The other parameters were taken to be the same 
as in the case of the steady state solution. The salient 
features of the time dependent solution are a transient at 
t = and a modulation of the intact fluorophore concen- 
tration at the driving frequency of the laser As expected 
the average level of intact fluorophores decreases with 
increasing average power, and the amplitude of the mod- 
ulation increases with increasing modulation amplitude. 
The approximate time dependent solution given by 
Eq. (36) was used to obtain the fluorescence signal given 
by Eq. (19). The fluorescence signal was detected by a 
computer model of a lock-in amplifier consisting of an 
in-phase and quadrature multipliers followed by low 
pass filters. The low pass filters were implemented by 
Fourier analyzing the signal, setting to zero all coeffi- 
cients corresponding to non zero frequencies, and then 
performing an inverse Fourier transform. The solid 
circles in Fig. 6a show the ratio of the detected quadra- 



ture and in-phase components at different modulation 
frequencies. Most of the parameters used in Eq. (36) 
were the same as in the steady state calculation except 
for width = 10 |um, F^ = 50 mW, and kj^ = 600 s"\ The 
solid line in Fig. 6a shows a fit to the calculated points of 
a fiinction given by R(f) = {a^f)l{b +/^ ), where/is the 
modulation frequency and a and b are fit parameters. 
Clearly the calculated response and R{f) form an excel- 
lent match. The solid points in Fig. 6b show experimen- 
tally measured values of the ratio of quadrature and in- 
phase response. The solid line in Fig 6b is a fit to the 
fiinction R{f). The match is reasonable. Thus there is a 
strong indication that the measured response can be rep- 
resented by the fiinction R{f), and that the fit parameters 
provide information about the various physical constants 
of the experiment. To elucidate the dependence of the fit 
parameters on the physical constants, it is best to use 
analytical solutions to Eq. (30) and obtain the analytical 
form of the response fiinction R{f). This task is carried 
out in a forthcoming paper. 



1.02 



1.00 - 




time, t 



Fig, 5, The time dependent numerical solutions of Eq. (30). Each curve corresponds to a different power of the laser beam. The ini- 
tial value is set to 1 in all cases. After an initial transient, the average concentration of intact fluorophores oscillates with the modula- 
tion frequency. The amplitude of the oscillation increases with increasing laser power. 
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Fig. 6. (a) The time dependent solutions, such as shown in Fig 5, were inserted into Eq. (19) to model the 
detected fluorescence signal. The fluorescence signal was inserted into a digital model of a two phase lock- 
in amplifier (discussed in the text) . The ratio of the outputs of the quadrature and in-phase lock-in amplifiers 
is shown by the solid circles. The solid line is a best fit to a function R(f) = {a^f)l{b +/^ ), where/ is the 
modulation frequency and a and b are parameters, (b) The sotid circles show the measured response of a flow- 
ing fluorescein solution hi a focused laser beam with a wavelength of 488 nm. The sotid line gives the best 
fit to the function R(f) given in (a). 
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We close this section with a brief description of the 
mathematical justification for the time dependent case. 
The reduction of the original system Eq. (4) and Eq. (5) 
to a single partial differential equation depended on the 
use of a multiple time scale perturbation expansion. 
The validity of the technique used for this problem 
(matched asymptotic expansions) is discussed in 
Appendix A. By spatially averaging the partial differ- 
ential equation we derive an ordinary differential equa- 
tion for the averaged reactant concentration as well as a 
simplified expression for the total fluorescence. These 
simplifications are possible because the varying func- 
tions (Eq. (21) and Eq. (24) respectively) are nearly 
constant. The proof of this fact can be found in 
Appendix B. Finally the equation for the spatial aver- 
age, as in Eq. (30), involves these functions and it is 
shown that replacing them by constants creates a small 
error. This is shown in Appendix C. 

6. Conclusion 

The dynamics of the populations of excited and 
unexcited fluorophores are described by a pair of first 
order partial differential equations that incorporate the 
kinetics of the fluorophore reaction driven by a period- 
ically varying laser beam and the convective fluid flow 
in a measurement apparatus. We obtain a dramatic 
reduction in the mathematical complexity of the model 
by taking advantage of the disparate time scales for 
excitation, relaxation and degradation. This is followed 
by averaging the reactant concentration over the laser 
beam, which is assumed to have a Gaussian profile. A 
single ordinary differential equation, Eq. (30), is 
obtained describing the evolution of the spatially aver- 
aged reactant concentration. The measured fluorescent 
signal, Eq. (19), is expressed in terms of the spatially 
averaged reactant concentrations. These equations pro- 
vide a prediction of the response whenever the photo- 
chemical reaction degrades the fluorophores flowing 
past the laser beam. 

Mathematical justification for the approximations 
used in the course of the derivation is given and some 
comparisons to the time independent case are 
described. The predicted frequency response (based on 
the time varying solution of Eq. (30)) is used to derive 
a functional form for the ratio of the in-phase and quad- 
rature parts of the fluorescent signal. This same form 
provides a good fit for experimental data. The fit and 
the connection with physical parameters of the experi- 
ment will be discussed in a forthcoming paper. 



7. APPENDIX A 

This appendix addresses the issue of the validity of 
the perturbation expansion of A^q ^^^ M for large times 
as stated in Eq. (6) and Eq. (7). The discussion is based 
on the standard observation that the solutions of Eq. (6) 
and Eq. (7) can be expressed in terms of the solution of 
a parameterized system of ordinary differential equa- 
tions. Here a new time variable s where t = s+0 is 
introduced. The relaxation rate is denoted by k^^. 

e^ie + s,y,e) = 

as 

£—-^(0 + s,y,£) = 

as 



K{x,y,e + s)N^{e + s,y,£) 
-k^N.iO + s,y,£)-£k,N,ie + s, y,£) 



(Al) 



dx 
ds 



(0 + S,£) = V 



N,iy,e,£) = N\ N,iye,£)=0, x(6>, £) ="4/. 



See for example the book Aris and Amundson [9], 
which contains many worked examples from the chem- 
istry literature as well as a general discussion of the 
transformation of partial differential equations such as 
Eq. (6) . This is known as the "method of characteris- 
tics." The variable 0>0 parameterizes the family of 
solution curves that satisfy the differential equations in 
s and the initial conditions at 5 = displayed on the last 
line. The original time variable is t= + s and this 
leads to the interpretation of as the time a particular 
family of particles entered the apparatus (at x = ~^/2) 
while s is the amount of time elapsed since entry. If we 
are following the set of particles that entered the appa- 
ratus at the beginning of the observation period, i.e., at 
0=0, then s = t. Here we will assume that > ~\ . For 
each 0, a perturbation method will be used to construct 
solutions of Eq. (Al) based on the representation, 

N,iy,e + s,£) = N^iy,e,s,£)+T],iye,T,£) 

N, (j, e + s,£) = N, (y, 0, s, £) +ri,(y 0, r, e) (A2) 

X(e + S,£)= X(e, S, £) + 8^(0, T, £) 
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where Ni / = 0, 1 refers to the large time component of 
TV, and rj , is the small time or inner correction of the 
large time solution. A similar decomposition holds for 
X. The inner correction governs the behavior of the 
solution of Al near the time of entry s = and since we 
expect a rapid initial transient it is written in terms of a 
fast time scale r = %. The large time components have 
perturbation expansions 

N,iy,e,s,e) = N,,iy,e,s)+eN,,iye, s) +0(8') (A3) 
X(0,s,e) = X^(0, s)+eX,(0, s)+0(£^). 



If we replace Nq,N^,x by their respective large time 
components in Eq. (Al) then, 



as 
as 



(A4) 



ds 



■ = V. 



On setting e = 0, one gets, 
Nioiy, 0, s) = = N^i ye, $ 



as 



(A5) 



the analogue of the quasi-equilibrium condition 
expressed in Eq. (9). Similarly when coefficients of £ 
on the left hand side of A4 are equated to those on the 
right, there results 



dN, 



ds 



ds 



oo_ _ 



= KN,,-K{X,{s),y,e+s)N,, 



^ dK{X,{s\y,e+s) ^^^ 
dX 



■■K{X,{s),y,e + s)N,, 



N„ 



-kK-kK + 



dx 



(A6) 



AJin 



ds 



1 _ 



= 0. 



Here we suppressed some of the arguments for ease 
of reading. Comparing this with Eq. (11) we see 
that for suitable initial conditions on Xq, N^ , and iVjo , 
N{y, 6, s) = iVoo (y> ^) ^) + Mo (f> ^) ^)) is a solution of 
Eq. (13) expressed in tenns of the solution of the para- 
meterized ordinary differential equation. 



dN{y,e,s) , k,kSX,{s),y,e + s) 



ds 



K+Kix^is\y,e+s) 

dX,{s) 



N(y,e,s) =0 



ds 



=v. 



The appropriate initial conditions can be obtained by 
examining the behavior of the inner corrections. 
Equations for them are obtained by returning to Al and 
replacing the functions Nq, Ni, and x and by the expres- 
sions on the right hand side of A2. We seek representa- 
tions of the corrections in terms of the following pertur- 
bation expansion: 

J7o (T, £, 6) = J]^ (T, 6) + £J7„i (T, 6) + 0(£') 

77, (T, £, e) = 7],, (T, e) + £77„ (T, 0) + 0(S^) (A7) 

^(r,e,e) = ^,(r,e) +e^,(r,e)+0(e') . 



We will suppress the dependence on y in what 
follows. After performing the previously described 
replacement involving Aland A2, and then using the 
right hand sides of A4, the desired equations are 
obtained. In them, the fast time variable is introduced 
by the change of variable, s = £T. We will write down 
the resulting equations just for rj^ and <^, 

^ = l(X(er)+e^(r),e+er)[N,(er)+no(r)] 

at 



di 

dr 



-k^(X(er) +e^(r), +er) N,(er) 
-UN,(eT)+ii,(T)]+k„N,(eT) 
-ek,[N,ier) +r],ir)] +ek,[N,ier) +r],ir)] 

= 0. 



(A8) 



The equation for tjq is quite similar to the one for rj^. 

The initial conditions of the original problem Al can 
be used as the initial conditions for rJQ^rj^, and ^ since 
these functions approximate the initial part of the solu- 
tion of Al. Thus ?7o(0)=A^', ^o(0) = 0, <^(0)="% 
Equations for the leading terms in the perturbation 
expansions in (A7) are then obtained by setting £ = in 
(A8). There results. 
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dri^ 



dz 
drj^ 

dx 
d^ 
dx 



= KiM^%0)ri,,iT)-k^ri,,iT) (A9) 



:0. 



The initial value Xq(0) is a constant that will be deter- 
mined shortly. Adding the first two equations in 

(A9) then gives us — (77^^ +7]^^) =0 so from the just 

dr 

stated initial conditions it follows that rj^J^r) +77io( ^) ~ 

r}QQ(0) + ri^Q(0) = N^ - that is, the sum of the leading 
terms in the expansion of the inner correction of 
A^o + A^i is a constant. The next step in determining the 
initial conditions for the large time solution (which in 
general will not be the initial conditions for the original 
problem) is to match the large time or outer solution 
with the inner correction in an overlapping region of 
the s axis. This region is defined as the outer reaches of 
the initial interval where the inner correction is valid 
and the beginning of the region of validity of the large 
time solution. To do this we calculate the limit of the 
inner correction as T ^ ^ by fixing s and letting £ ^ 0. 
The limit must match the limit of the large time solu- 
tion as s = er^O where t is fixed and £ ^ 0. For our 
problem this is written as, 

N' = lim(r7„„(T) +r7,„(T)) =\im( N^is) + N,,is)) 



-^ = lim<^„(T)=lim X,(s). 



(AlO) 



In this way we determine the initial condition for the 
leading order large time solution N(y,0,s) and thereby 
determine the boundary condition. Indeed from the 
second equation in AlO, the initial condition for the x 
position to leading order is Xq(0) = "%. Given the large 
time solution and inner corrections formally construct- 
ed here, we can for each 0>0 construct a solution of 
Al that is uniformly valid throughout the interval 
< s < ^/£ for each fixed £ > Standard theorems in 
the theory of singular perturbations show that the oper- 
ations used here to obtain it are rigorously justified. The 
key hypothesis that must be satisfied in our problem is 
that the inner solution must approach a finite limit 
exponentially. Recall that ti^q = N^- tjqq. Substituting 
this relation into the first equation in A9 results in the 
following equation for r]oo- 

^HKi^oiO%0)+kJri^=k^N' (All) 



dr 



Solving this constant coefficient equation shows that 
r]oo and thus r]io has the required behavior as T— > ^. 
Note that equation (Al) is not a standard singular per- 
turbation problem because the Jacobian associated with 
the reduced problem, i.e., the problem obtained from 
setting £ = in the first two equations of (Al), is singu- 
lar. Validity of the singular perturbation method in this 
"critical" case was proved by Vasil'eva and Butuzov 
[10] and was treated analytically and numerically by 
O'Malley and Flaherty [11]. These results require that 
the Jacobian have a constant deficient rank k for all s 
and its non-zero eigenvalues have negative real parts. A 
discussion of these problems can be found in O'Malley 
[12] and Smith [13]. 

A complete solution of (Al) uniformly valid in an 
interval comes from combining the large time and inner 
correction and then subtracting the common value they 
share in the overlapping region. lfN(y,0,s) is the total 
fluorophore population to leading order we have. 



N^(y,e + s,£) = ^ 



k^Niy^O^s) 



kSX,{e,sly,e + s) + k^ 



+ N'e 



k+k 



-(*«+*« )y 



-A^" 



kN" 



(A12) 



- = " - +0{e) 
k+k„ 



where k^ =k^(X,(0),y,e) =k^(-y^, y,e). Note that 
when s is large 



N,iy,0 + s,e) = - ^ J^^^^'^'^l - + 0(e) + d(s) . 



kSX,(s),y,e + s) + k, 
Similar reasoning shows that 



(A13) 



^oo(0)=^ 



- mo 



k„+k^iX,is),y,0 + s) (A14) 
+ 0(s) + l}(s) 

where 5 (s) and P (s) are exponentially small in s and 
0(e) holds uniformly in the interval 0<s <^/£ and 
where Xq(s) = ~^/2-^ vs. On adding Eqs. (A13) and 
(A 14) we see that at large times, the total fluorophore 
population is dominated by N(y,0,s). It is in this sense 
that we can justify the statements in Sec. 2. 
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8. APPENDIX B 



We will write N(r,t) = 0(r) +^(r, /) where 0( r) 

1 ^ 
lim — {N(r, t) (i/and ^(r, /) = N{r, t) - 0( r). We begin 



will have passed and we have continuity in ^. The 
solution of Eq. (B2) for t > ^^is 



' 1 r [ 11- 1 

byshowingthat^(n0and3^(^'^){^ are small when ^(r,0=- ^xp|--J//(^,;.,T>/^ |g (x>,t)^^ 



the parameter r\ is small. To obtain an expression 
for ^ we will first derive an equation for O from 

1 ^ 
Eq. (20), by applying the operator lim— { }. 

One finds O that must satisfy ^ 



aO(r) P f(r) 



(Bl) 



. L/ < y < L/ 

72-^-/2 
0(r)| 



-a<y<a 



,= N' 



We assume that all of the limiting averaging opera- 
tions used to derive B 1 are valid. This will be the case 
for example when the integrands are quasi-periodic 
[14] - a reasonable assumption for this problem. The 
boundary condition at x = I2, comes from applying the 

1 ^ 
operator lim — U } to the boundary condition Eq. (15). 

Substituting Eq. (Bl) into Eq. (20) will lead to an 
equation for *?. 

^ '+v y ' +H(r,t)^'{r,t)=g{r,t) 



dt 
H{r,t) = 



dx 
\ + bP{r,ty 



g{r,t) = 



( 



^(r,OI 



-iHr,t)^ 



r\PJir) 
^ + bPJ{r) 



\ 



0(r) 



-Y2 



= , ^(?;0) =-0(r) , x>-L/, 



(B2) 



The conditions for ^ at x = 1 2, are obtained from 
the conditions for O and A^. Note the discontinuity in 
the value of ^ at / = as particles are pumped into 
the vessel. It propagates to larger values of x for / > as 
a shock. However since here we will assume that 



t>- 



(x+A/) 



— for any x (for example if / > V the shock 



(B3) 



x-^ 



where r = (x,^, /) = /-I - — — 1. The derivative with 

I ^ ) 
respect to x is written, 

dx 

2 

+g{r,t) - ^ J exp -\\HdB, \g,dx' 

T 

(B4) 

where for easier reading we have suppressed some of 
the function arguments. The functions H^ and g^ are 
derivatives 

From Eq. (B2) it can be seen that g^ = -H^ \^=^^ 0(r). 



To show that | ^(r, /) | and 



dx 



are small we will 



show that g, //, H^ are all small. Indeed from Eq. (Bl) 
we have 0(r)<A^^ so |gj will be small if |//^| is. 
Using the fact that P(rj) = (Pq + A(t)) f (r) where 
A(/) = P^cos(o)t) we find that 

dHjrj) ^ rif(r)Q)P,sm(m) ^^^^ 
dr [l^bP(rj)f(r)]' ' 

We can perform a similar analysis of the remaining 
terms in Eq. (B3) and Eq. (B4). Using the notation 0(p) 
to denote a term whose absolute value remains bound- 
ed by a multiple of p as p ^ 0, the results are summa- 
rized as 
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3//(r,T) 
3t 

3t 
//(r,T) =0(?]) 



0(?]) 



0(?]) 



(B7) 



??A(0/(r) 



(l + Z./;/(r))(l+^P(r,0/(r)) 



0(?]) 



These relations hold uniformly for all ^ > since all 
the time varying functions in (B7) are periodic. Thus by 
considering 7] small enough we can conclude that for 
alU>0 



^(r,t) = 0(r]X 



dx 



= 0(ri) 



Nir,0 = Oir).O(,), milAJ.m^oi,)}'''^ 



dx 



dx 



Here we used the relation N(rj) = <|)(r) + ^(rj). 

Finally we turn to a discussion of the replacement of 
the functions a(t) and aXt) in Eq. (22), Eq. (24), and 
Eq. (25) by constants a and a\ We use the fact that the 
constants are 



a 



R ^i dx 



(B9) 



to write a(t)= a+ a (t) where 



JJ/o \ijm\- ijmlijn' 



«(0 = - 



jjf[o+w]\jjfo 



■ (BIO) 



Applying Eq. (B8) for ^ to the numerator of (BIO) we 
see that 



a(t) = a-\-0(ri). 



In an analogous fashion the function a'Q) = a'+aXt) 
where 



aXt) = 



R 9-^ Kr 9-^ 



dx 



iln 



dx dx •'/ dx 



(B12) 



3^ 

Using Eq.(B8) for and applying it to Eq.(B12) 

dx 

therefore implies the equation 

a\t) = a+0(ri). (B13) 

Equations (Bll) and (B13) hold uniformly for all ^ > 0. 
Thus the effect of replacing a(t) by the constant a is 
to produce an error of size 0(7]) by Eq. (Bll) and 
Eq. (B13). This argument also justifies the replacement 
in Eq. (19). 

9. APPENDIX C 

In Sec. 3 an equation Eq. (26) for the evolution of the 
spatially averaged reactant concentration is derived 
from Eq. (20). To properly assess the effect of replac- 
ing the functions a(t) and aXt) respectively with con- 
stants a and a' in the derivation of Eq. (26) it must first 
be recognized that the time derivative of N(r,t) is small 
since it is 0(7]). This can be seen from the definition of 
^ in Appendix B and Eqs. (B2), (B7), and (B8). We 
must therefore show that the neglected terms are of 
smaller asymptotic order. Before the replacement, the 
left hand side of Eq. (20) after multiplication hyf(x,y) 
and integration is actually Eq. (22). Using the notations 
of Appendix B, Eq. (22) is rewritten, 



d{m)) 



dt 



+ b(P,+A(t)) 



■ [l+b(P,+A(t))a] 



|{(«.a(0)j(J 



f(r)N(r,t)dr\ 



dt 



.KP„^A(0)j[^](MO)-(0^}. (CI) 



The term a(t) 



d{N(t))\ 



dt 



is O(J70 thanks to Eqs. (B8) 



(Bll) b(P,+A(t)) 



and (Bll) and can therefore be neglected. The term 

3a(0 



dt 



{N(t)) is proportional to 7]bI^(sGG 



Eq. (B5), Eq. (B6), and Eq. (BIO)). Thus this term can 
be dropped when either P| or b is small. Note that if the 
constant ais large, (i.e., a» 1) the following terms are 
retained because they are going to be of larger order 
than the terms that are dropped. 



b(P,-^A(t))a 



d{N(t))]_ 



dt 



= 0(7]bP,a) 



-7]a(P,+A(t)){N(t)) = 0(7]a) 
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The first expression comes from the left hand side of 
C 1 and it is retained, and the second expression comes 
from the right hand side of Eq. (20) after multiplication 
and integration, i.e., Eq. (25). Therefore let us write 
the right hand side of Eq. (20) after multiplication by 
fix,y) and integration as. 



dr 



(C2) 



-v(l + biP, + A(0)(a' + a'(0)« JJ /( r) ^^^ 

-ri(a-^a(t))(P,+A(t)){N(t)). 

We observe that 

-r7a«(/'o+A(0)(A^(0) = 0(r7^) 



because of Eq. (Bll) and Eq. (B13). Again taking b small 
means that both these terms can be dropped. Moreover 

if a'» 1, the term biP^ +A(0) a'\\ f(r) M^!!^ dr is 

JJ dx 

R 

retained. Thus, the remaining terms in Eq. (CI) form 
the left hand side of Eq. (26) while the remaining terms 
in Eq. (C2) become the right hand side of Eq. (26). 
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